import struct
import bz2
import numpy as np
from matplotlib import pyplot as plt
import matplotlib


def show_hot(rs):
    # rs = np.maximum(rs, -100)
    # rs = np.minimum(rs, 2560)
    # cmap = matplotlib.cm.jet
    # # 设定每个图的colormap和colorbar所表示范围是一样的，即归一化
    # norm = matplotlib.colors.Normalize(vmin=0, vmax=950)
    # # 显示图形，此处没有使用contourf #>>>ctf=plt.contourf(grid_x,grid_y,grid_z)
    # gci = plt.imshow(rs, cmap=cmap, norm=norm)
    # plt.colorbar()
    # plt.show()

    rs = np.clip(rs, -3000, None)
    plt.imshow(rs)
    plt.colorbar()
    plt.show()


def readFile(filename):
    # print("read file :" + filename)

    fd = open(filename, 'rb')
    headInfo = {}
    headInfo['flag'] = fd.read(4).decode('utf-8')  # 0-3	文件固定标识：MOC
    headInfo['version'] = fd.read(4)  # 4-7	文件格式版本代码，如: 1.0，1.1，etc
    headInfo['intfileBytes'] = struct.unpack('i', fd.read(4))[0]  # 8-11	包含头信息在内的文件字节数，不超过2M
    headInfo['productCode'] = struct.unpack('h', fd.read(2))[0]  # 12-13 拼图产品编号
    headInfo['coordinateSystem'] = struct.unpack('h', fd.read(2))[0]  # 14-15 坐标类型: 2=笛卡儿坐标,3=等经纬网格坐标
    headInfo['productName'] = fd.read(8)  # 16-23 产品代码,如: ET,VIL,CR,CAP,OHP,OHPC
    headInfo['productDesc'] = fd.read(64)  # 24-87 产品描述,如Composite Reflectivity mosaic
    headInfo['intblockPositon'] = struct.unpack('i', fd.read(4))[0]  # 88-91 产品数据起始位置(字节顺序)
    headInfo['intblockBytes'] = struct.unpack('i', fd.read(4))[0]  # 92-95  产品数据字节数

    headInfo['inttimeZone'] = struct.unpack('i', fd.read(4))[0]  # 96-99 数据时钟,0=世界时,28800=北京时
    headInfo['year'] = struct.unpack('h', fd.read(2))[0]  # 100-101	观测时间中的年份
    headInfo['month'] = struct.unpack('h', fd.read(2))[0]  # 102-103	观测时间中的月份（1－12）
    headInfo['day'] = struct.unpack('h', fd.read(2))[0]  # 104-105	观测时间中的日期（1－31）
    headInfo['hour'] = struct.unpack('h', fd.read(2))[0]  # 106-107	观测时间中的小时（00－23）
    headInfo['minite'] = struct.unpack('h', fd.read(2))[0]  # 108-109	观测时间中的分（00－59）
    headInfo['second'] = struct.unpack('h', fd.read(2))[0]  # 110-111	观测时间中的秒（00－59）
    headInfo['intobserveSecond'] = struct.unpack('i', fd.read(4))[0]  # 112-115	观测时间的seconds[年月日分秒,转换成秒的单位]
    headInfo['observeDate'] = struct.unpack('H', fd.read(2))[0]  # 116-117 观测时间中的Julian dates
    headInfo['productDate'] = struct.unpack('H', fd.read(2))[0]  # 118-119 产品处理时间的天数[花了几天时间解析出来]
    headInfo['intproductSecond'] = struct.unpack('i', fd.read(4))[0]  # 120-123 产品处理时间的描述

    headInfo['intymin'] = struct.unpack('i', fd.read(4))[0]  # 124-127 数据区的南边界，单位：1/10000度，放大1万倍
    headInfo['intxmin'] = struct.unpack('i', fd.read(4))[0]  # 128-131 数据区的西边界，单位：1/10000度，放大1万倍
    headInfo['intymax'] = struct.unpack('i', fd.read(4))[0]  # 132-135 数据区的北边界，单位：1/10000度，放大1万倍
    headInfo['intxmax'] = struct.unpack('i', fd.read(4))[0]  # 136-139 数据区的东边界，单位：1/10000度，放大1万倍
    headInfo['intxCenter'] = struct.unpack('i', fd.read(4))[0]  # 140-143 数据区中心坐标，单位：1/10000度，放大1万倍
    headInfo['intyCenter'] = struct.unpack('i', fd.read(4))[0]  # 144-147 数据区中心坐标，单位：1/10000度，放大1万倍
    headInfo['intxCount'] = struct.unpack('i', fd.read(4))[0]  # 148-151 格点坐标为列数,
    headInfo['intyCount'] = struct.unpack('i', fd.read(4))[0]  # 152-155 格点坐标为行数,
    headInfo['intxGap'] = struct.unpack('i', fd.read(4))[0]  # 156-159 格点坐标为列分辨率，单位：1/10000度，放大1万倍
    headInfo['intyGap'] = struct.unpack('i', fd.read(4))[0]  # 160-163 格点坐标为行分辨率，单位：1/10000度，放大1万倍
    headInfo['height'] = struct.unpack('h', fd.read(2))[0]  # 164-165 雷达高度[距离地面距离]
    headInfo['compress'] = struct.unpack('h', fd.read(2))[0]  # 166-167 数据压缩标识, 0=无,1=bz2,2=zip,3=lzw[压缩分类]
    headInfo['intradarCount'] = struct.unpack('i', fd.read(4))[0]  # 168-171有多少个雷达进行了拼图
    headInfo['intunzipBytes'] = struct.unpack('i', fd.read(4))[0]  # 数据段压缩前的字节数
    headInfo['scale'] = struct.unpack('h', fd.read(2))[0]  # 176-177
    headInfo['shortunUsed'] = struct.unpack('h', fd.read(2))[0]  # 178-179
    headInfo['RgnID'] = fd.read(8)  # 180-187
    headInfo['units'] = fd.read(8)  # 180-187
    headInfo['reserved'] = fd.read(60)  # 188-255留用

    # print(headInfo)

    bi = fd.read()
    data = bz2.decompress(bi)
    ds = struct.unpack('26040000h', data)
    array = np.array(ds, dtype=np.int16)
    rs = np.reshape(array, (4200, 6200))

    return {"headInfo": headInfo, "data": rs}


def lon_lat_bounds(data):
    lon = (data['headInfo']['intxmin'] / 1000.0, data['headInfo']['intxmax'] / 1000.0)  # 经度longitude
    lat = (data['headInfo']['intymin'] / 1000.0, data['headInfo']['intymax'] / 1000.0)  # 纬度latitude

    print("经度lon:", lon, "纬度lat:", lat)

    return (lon, lat)


def normalization(data):
    data = np.clip(data, 0, None)
    data = data.astype(np.uint16)
    return data


if __name__ == '__main__':
    # files = os.listdir('CREF')
    #
    # i = 0
    # for file in files:
    #     filename = "CREF\\" + file
    #     data = readFile(filename)
    #     # plt.imshow(data)
    #     # if i == 0:
    #     #     plt.colorbar()
    #     # i += 1
    #     # plt.savefig("out\\" + file + '.png')  # 先存，再show
    #     # plt.show()

    # file = "CREF/ACHN_CREF_20171113_000000.bin"
    # data = readFile(file)
    # print(data['headInfo'])
    # # lon_lat_bounds(data)
    # show_hot(data['data'])

    import config
    radarInfo = readFile(config.CREF_FILES + 'ACHN_CREF_20171113_035000.bin')  # 读取雷达文件
    lon_lat_bounds(radarInfo)
    # show_hot(rd)
